Multiconfigurational Hartree-Fock theory for identical bosons in a double well 
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Multiconfigurational Hartree-Fock theory is presented and implemented in an investigation of 
the fragmentation of a Bose-Einstein condensate made of identical bosonic atoms in a double well 
potential at zero temperature. The approach builds in the effects of the condensate mean field and 
of atomic correlations by describing generalized many-body states that are composed of multiple 
configurations which incorporate atomic interactions. Nonlinear and linear optimization is utilized 
in conjunction with the variational and Hylleraas-Undheim theorems to find the optimal ground and 
excited states of the interacting system. The resulting energy spectrum and associated eigenstates 
are presented as a function of double well barrier height. Delocalized and localized single config- 
urational states are found in the extreme limits of the simple and fragmented condensate ground 
states, while multiconfigurational states and macroscopic quantum superposition states are revealed 
throughout the full extent of barrier heights. Comparison is made to existing theories that either 
neglect mean field or correlation effects and it is found that contributions from both interactions are 
essential in order to obtain a robust microscopic understanding of the condensate's atomic structure 
throughout the fragmentation process. 

PACS numbers: 



I. INTRODUCTION 

Recent experimental realization of a trapped atom in- 
terferometer using a coherently split Bose-Einstein con- 
densate (BEC) 12 as well as the direct observation of 
tunneling and self-trapping in weakly linked BECs H H 
have provided an impetus to formulate a comprehensive 
theoretical description of the zero temperature Bose gas 
confined to a trapping potential that can be continuously 
deformed from a single well into a double well with large 
barrier height. Such a deformation of the trapping po- 
tential impels the BEC from a single coherent entity to a 
fragmented condensate made up of two coherent and po- 
tentially correlated moieties. The quantum many-body 
physics governing this complex fragmentation process in- 
volves competition and balance between the effects of the 
condensate's mean field on the interacting atomic gas and 
the correlations that emerge between atoms in different 
Fock states. 

Various theoretical descriptions of the simple and frag- 
mented BEC exist in the literature today. At both low 
and high barrier heights, the many-body ground state of 
the condensate is well approximated in mean field the- 
ory by a single Fock state, which expresses a particular 
arrangement of atoms among one, two, or even many 
single- particle states 0, IE ■ ^ OIU y one single-particle 
state is involved, then mean field theory reduces to the 
Hartree theory [jj, which further reduces to the stan- 
dard Gross-Pitaevskii formalism l9Ul0ll upon invoking the 
contact interaction approximation [llL . More elabo- 
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rate mean field theories, such as Hartree-Fock [l3|], are 
built by utilizing more than one single-particle state and 
imposing the correct symmetrization due to the indis- 
tinguishability of identical particles. Where mean field 
equations specify the underlying single-particle states, we 
call the single Fock state a single configuration . How- 
ever, more complex many-body states exist, at all barrier 
heights, that are made up of superpositions of multiple 
configurations. 

As mean field theory describes only a single config- 
uration, it necessarily lacks all correlation effects that 
arise between configurations. To this end, multiconfig- 
urational approaches have been attempted, which when 
applied to cold atomic gases in a double well trap ge- 
ometry, are largely based on two-well limits of contin- 
uum lattice models such as the Bose-Hubbard model. 
Both atomic correlations and, to a partial extent, mean 
field interactions are included automatically in these the- 
ories. Correlations emerge between configurations while 
mean field interactions occur in two places: first, di- 
rectly through the interaction term in the many-body 
Hamiltonian, and second in the underlying single-particle 
states that make up the matrix elements in the Hamil- 
tonian as well as each configuration. Previous multi- 
configurational efforts found in the BEC literature in- 
clude only the first type of mean field effect and there- 
fore do not describe general many-body states of the sys- 
tem by superpositions of configurations but only those 
where the effects of atom-atom interactions on the shape 
of the single-particle states are neglected. In these works, 
the underlying single-particle states have been chosen 
either to be solutions of the single-particle Schrodinger 
equation |l4j or parameters have been introduced to re- 
place the matrix elements in the Hamiltonian altogether 
111 El O E El For example, in the Bose-Hubbard 
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model U is a site energy and J is a tunneling param- 
eter, and neither depends on N, the particle number, 
and neither is computed from first principles. In other 
words, mean field effects arc not fully included in these 
treatments: the single-configurational and multiple Fock 
state approaches are separate and complementary. 

A more complete description should characterize a gen- 
eral state of the system by a superposition of many con- 
figurations in which the underlying single-particle states 
include the effects of the condensate mean field. Such 
a theory would be capable of illuminating the quantum 
many-body structure of the BEC throughout the frag- 
mentation process as well as at the extreme simple and 
fragmented limits. It is the purpose of the present paper 
to provide such a formulation. 

The challenge to formulate such a theory has alre ady 
been partially fulfilled by several authors 0, IS LL-I OH 
EiJ E3 EB E3 However, to our knowledge, no satisfac- 
tory work has been presented in the cold atom context 
that fully addresses this task. To this end, we formulate 
a new approach that variationally combines the Hartree- 
Fock mean field theory for N identical bosons in two 
single-particle states with a full diagonalization of the 
many-body Hamiltonian restricted to a basis of N + 1 
generalized configurations stemming from each Hartree- 
Fock configuration. This allows for a description of the 
Bose gas at zero temperature where atomic correlations 
emerge between configurations into which mean field ef- 
fects are built. Due to its composition, this approach 
incorporates both types of interaction and serves as a di- 
dactic device for elucidating where they become impor- 
tant in the fragmentation process as well as how these 
two intertwined but distinct effects change the system. 
By choosing the atomic interaction strength to be zero 
in the bosonic Hartree-Fock mean field equations, our ap- 
proach recovers the Schrodinger based model developed 
in |l4j . while limiting our generalized many-body state 
to a single configuration recovers the mean field theories 
of 0, U A preliminary investigation of our formalism 
can be found in |2fl |. We acknowledge that a fermionic 
analog of our model, called multiconfigurational self- 
consistent field theory, is widely known in quantum chem- 
istry and has been quite successful in accurately describ- 
ing atomic and molecular electronic structure andgeom- 
etry both at equilibrium and at dissociation 0, |22, H^J • 
In fact, a proper description of dissociation of poly- 
atomic molecules in the Born-Oppenheimer approxima- 
tion [2lL I22I l23l | is closely related to the problem of frag- 
mentation of a BEC into two or more fragments. Related 
time-independent 0> E3 an d time-dependent j2|| mul- 
ticonfigurational approaches have also been developed to 
treat molecular vibrations at the Hartree level. In the 
spirit of these efforts, we refer to the work developed in 
this paper as the multiconfigurational bosonic Hartree- 
Fock theory. 

To establish a consistent and general enough nota- 
tion, which can be extended to our multiconfigurational 
Hartree-Fock approach, we organize the paper as fol- 



lows. In Section II, we review the many-body theory 
of a gas of N identical bosons restricted to a finite model 
space. Ground and excited eigenstates of the system are 
expanded onto a basis consisting of N + 1 Fock states 
made up of two single-particle states. Model calcula- 
tions, which lack the effects of atomic interactions on the 
underlying single-particle states, are discussed where the 
trapping potential is deformed from an initially single 
well to a double well geometry. Section III is devoted to 
a survey of Hartree-Fock mean field theory for N identi- 
cal bosons in two single-particle states. It is emphasized 
that mean field theory describes only a single configura- 
tional state and therefore lacks the correlations described 
by superpositions of multiple configurations. Imaginary 
time integration is briefly discussed for efficient solution 
of the resulting coupled nonlinear differential equations. 
While imaginary time integration schemes are a stan- 
dard method of solution for the bosonic mean field equa- 
tions with one single-particle state, i.e., for the Gross- 
Pitaevskii equation 27], or with multicomponent spinors 
|28| , we highlight our method of maintaining spatial or- 
thogonality between two single-particle states. Solutions 
of the mean field and Schrodinger equations are com- 
pared at various barrier heights in the strongly interact- 
ing limit. 

In the Section IV, we introduce the multiconfigura- 
tional bosonic Hartree-Fock theory. For each individual 
energy level, multiconfigurational Hartree-Fock states of 
the system are constructed from the variationally opti- 
mal linear combination of generalized configurations in 
which the underlying mean field states are chosen so that 
the partitioning of atoms between its two states allows 
the corresponding energy to be minimized. Lastly in Sec- 
tion V, a systematic investigation of the atomic structure 
of the BEC, as a function of barrier height, throughout 
the fragmentation process is carried out. This is fol- 
lowed in Section VI by a summary and an indication 
of further work needed to describe ongoing experiments. 
An Appendix is devoted to an informal statement of the 
Hylleraas-Undheim theorem, which justifies such a state- 
by-state use of the variational theorem. In particular, 
we discuss our application of this theorem to the opti- 
mization of BEC excited states and their corresponding 
energies. 



II. REVIEW OF MANY-BOSON THEORY IN A 
RESTRICTED BASIS 

The many-body Hamiltonian for a gas of N identical 
spinless bosonic atoms of mass m at zero temperature 
interacting via a two-body potential V^(x, x') = V(\x — 
x'|) is given by (2i| 

H = /*t( x ) {(_n 2 /2m)V 2 + VWt(x)} #(x)d 3 x 

+ (l/2)/* t (x)# t (x')^(x,x')4'(x')*(x)d 3 xd 3 x'. 

(1) 
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Here iff and iff' are boson field operators which satisfy 
the equal time commutation relations [V(x), ^(x')] = 
5(x-x') and [*(x),*(x')] = [* f (x), £t( x ')] = q. The 
external trapping potential V ext that we have in mind 
throughout this paper is one that can be continuously 
deformed from a single well to a double well geometry. 

A. Restriction to Two Single- Particle States 

The atomic structure of a BEC confined in a double 
well trapping potential at zero temperature can be rea- 
sonably described in a basis of restricted Fock states of 
the form BUI 

\N U N 2 ) = ®)^)^vac)/VWW, (2) 

where Nk atoms are in each of the single-particle states 
\Xk) = bj.\vsLc) for k = 1,2 and the total number of 
atoms is fixed at N — N\ + N 2 . Here, |vac) is the vac- 
uum state in which no atoms are present. The operators 
St and bk are boson creation and annihilation operators 
that add and remove single atoms in the \xk}- They sat- 
isfy the basic commutation relations [6fe,St] = Ski and 
[b k ,bi] = [b{,b\] = for k,l = 1,2. While © is cer- 
tainly not the most general eigenstate imaginable, we 
are interested mainly in the zero temperature proper- 
ties of the condensate as its constituent atoms are ex- 
changed between two macroscopically occupied single- 
particle states. Where the effects of finite temperature 
and of fragmentation into more than two condensates are 



negligible, the states provide a rich basis in which to 
explore the many-body physics of the BEC fragmentation 
process. 

The set of all Fock states of the form (J2J, with all pos- 
sible numbers of atoms in each of the two single-particle 
states, exhausts the restricted iV-boson Fock space. That 
is, the model space is spanned by the collection 

{|AT,0>, lJV-1,1), |iV- 2,2), \0,N)}, (3) 

which is taken as a complete set having N + 1 elements. 
Therefore, it is possible to expand an eigenstate of the 
many-body Hamiltonian Q as a linear combination of 
N + 1 individual Fock states according to [lj, [lj| 

JV 

|*%= C* 1 \N 1 ,N 2 = N-N 1 ), (4) 

ATi=0 

where the expansion coefficient C%- expresses the prob- 
ability amplitude for the ^th excited state of the system 
to be in |iVi, JV 2 ). 

The many-body Hamiltonian that is associated with 
this model space may be derived from Q by substitution 
of the two-state expansion of the boson field operator 

*(x) = X i(x)6i+X 2 (x)6 2 , (5) 
where the expansion coefficients Xfe( x ) = ( x |xfc) are CCH 
ordinate space single-particle wavefunctions that are, as 
yet, unspecified. After some standard algebra, one ob- 
tains m 



H = huNi + h 22 N 2 + (l/2)[Fi2i2 + Vi 2 2i]AiA 2 + (1/2)[V 2121 + V 2112 ] N 2 N X + [h 12 + Vm 2 (iVi - 1) + V 2221 N 2 ]b\b 2 

+ [h 2 i + V xll2 N x + V 2221 (N 2 - l)]blh + {l/2)[V U u(Nl - N ± ) + V 2222 {N% - N 2 ) + V 1122 (b\b\b 2 b 2 + & 2 & 2 &i&i)] 

(6) 

I 



where Nk = b k bk is the occupation number operator for 
each state |xfc), N = N x + N 2 is the total particle number 
operator, and 

hm = J Xfe(x){(-n 2 /2m)V 2 + V cxt (x)} X i(x)d 3 x 

Vkimn = /Xfc(x)x/(x')V A (x,x')x m (x)x„(x')d 3 xd 3 x' 

with k,l,m,n = 1,2 are matrix elements of the single- 
particle Hamiltonian /i(x) and two-body interaction po- 
tential V(x, x'). Throughout this paper, the wavefunc- 
tions Xi( x ) and X2( x ) are real- valued functions so that 

Vklrnn Hfcnm Vrnnkl Vnrnlk 'mlkn 
— Vlmnk — Vknml — Kifcim- 

This Hamiltonian is more general than that, for ex- 
ample, of Spekkens and Sipe [T^j in that the external 



potential T4 x t is not assumed to be symmetric a priori 
since we anticipate the possibility for deformation of a 
single well into either a symmetric or asymmetric dou- 
ble well. Exclusion of certain nonlinear terms and the 
assumption of N- independent hki and Vkimn , reduces © 
to a two-state Bose-Hubbard Hamiltonian. In the ther- 
modynamic limit, Bose-Hubbard theory provides, inter 
alia, the standard model for the description of zero tem- 
perature quantum phase transitions |3fll l3l| in atomic 
gases confined to optical lattices |32|. 

With the Hamiltonian © and associated eigenstates 
0, the many-body Schrodinger equation 

H\* N ) V = E»\* N ) V , (9) 

which is a linear equation, may be represented in the 
restricted Fock state basis ©. This results in the set of 
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(N + 1) x (N + 1) matrix equations 
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E v (10) 



with matrix elements and overlap 

H IilNl = {N[,N-N' 1 \H\N l ,N-N 1 ) 
8 N[Nl = {N[,N — N[\Ni,N — N%) 

for N[,Ni — 0, . . . , N. Once the wavefunctions %i and 
Xi are specified, the matrix elements hki and Vkimn are 
defined and can be used to build Hn'Ni ■ Then, the eigen- 
value equations (fTUfl can, in principle, be solved to give 
the expansion coefficients C Ni and energies E v of the 
ground state and first TV excited states of the system. 
This computation may be repeated as the external po- 
tential V^ x t is deformed from a single well to a double 
well where the geometry encourages break up into two 
condensate fragments. Following Penrose and Onsager 
|33|. a fragmented condensate may be characterized by 
the presence of at least two large eigenvalues of the single- 
particle reduced density matrix 0113 

7 "(x,x') = ^^(x^x')!^. (12) 

Within this model space, the full diagonalization of 
(|1U|I has been found, in practice, to be achievable for sys- 
tems with a large number of atoms. This is facilitated 
by the significant reduction of the Hamiltonian matrix in 
HI 0| l to a pentadiagonal form . Standard linear alge- 
bra routines j3|J are well suited for the resulting banded 
eigenvalue problem. 

B. Model Calculations 

The energies of the ground and first N excited states of 
the Hamiltonian 10 have been mapped out as a function 
of the single dimensionless parameter a in correlation 
diagrams by Reinhardt and Perry as illustrated in 
Figure □ In [3, it has been assumed that the single- 
particle energies and matrix elements of the two-body 
potential V, in the standard contact interaction approx- 
imation V(x, x') = gS(x — x'), can be parameterized ac- 
cording to [53 

Vkkkk — hkk — fuJ 

Vkkki = hki = -^exp(-a) (13) 
Vkku = hajexp(-a) 

in terms of the harmonic oscillator energy Huj and length 
f3 = y/h/mu) for a symmetric double well trapping po- 
tential with k 7^ I = 1,2. Variation of the parameter a 
allows for a continuous change between strong tunneling 
(small a) and weak tunneling (large a) regimes. Within 




FIG. 1: (Color online) Parameterized model ground and ex- 
cited state energies for g = 0.1 hw ■ /3 3 and iV = 20 as a func- 
tion of a. Note the energy level mergings and resulting ridge 
structure separating BEC (small a) and fragmented states 
(large a). The Fock states below the ridge are delocalized 
and nondegenerate while the states above the ridge are local- 
ized and doubly degenerate. The alternation of line style has 
been chosen to aid in visualization. 

this simple ansatz, where none of the parameters depend 
upon the particle number N, the wavefunctions xi and 
Xi are not specified. The ridge structure in Figure ^ 
shows the boundary between simple BEC for small a 
and fragmented BEC for large a. When working in the 
Fock basis © where N\ and N2 are the number of atoms 
localized in the left and right wells of a double well po- 
tential, it has been found in |l9j | that the distribution of 
Fock states contributing to the ground state below the 
ridge is binomial in form, while the distribution of Fock 
states above the ridge reveals the existence of macro- 
scopic quantum superposition states. The study of the 
correlation diagrams in Figure ^ has led to the predic- 
tion that many interesting highly excited states, such as 
Schrodinger cats, exist in the weak tunneling regime and 
may be created by phase engineering [l^, |33 • 

Other authors [14J opt to approximate the single- 
particle wavefunctions xi and xi f° r a given barrier 
height by solutions of the single-particle Schrodinger 
equation 

Mx)Xfc(x) =£fcXfc(x) (14) 

where £k = hkk/ (Xk\Xk) is the single-particle Schrodinger 
energy for state k and where the wavefunctions are in- 
dependent of the particle number N and value of g. The 
two-state Hamiltonian JBJ, which does depend upon N 
and upon g in the contact interaction approximation, is 
then diagonalized. Figure [21 presents the associated en- 
ergy eigenvalues as a function of barrier height, where 
Xi and xi are computed from (|14|) . A qualitatively sim- 
ilar energy correlation diagram exists with ridge struc- 
ture marking the phase transition between nondegener- 
ate states and doubly degenerate states. This approach is 
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barrier height (Hui) 

FIG. 2: (Color online) Schrodinger model ground and excited 
state energies for g — 0.1 huj ■ (3 3 and N — 20 as a function of 
barrier height. Note the qualitative similarities of this energy 
level correlation diagram to that of Figure The solid curve 
beginning at 20 hw is the single-particle Schrodinger energy 
while the remaining curves are the eigenvalues of the Hamil- 
tonian © built from the Schrodinger solutions at each barrier 
height. Note that the Schrodinger energy iViei + iV2£2 differs 
from the many-body eigenvalues E" because it includes no 
atom-atom interaction. The alternation of line style has been 
chosen to aid in visualization. 



justifiable for weakly interacting atomic gases where the 
effect of atomic interactions on the shape of the wave- 
functions is small, however, it breaks down wherever 
atom-atom interactions are important enough to affect 
the value of the parameters or matrix elements them- 
selves. This is the case in the experiments discussed in 
[1I1H0, where Thomas-Fermi mean field effects domi- 
nate the shape of the atomic wavefunctions in a strongly 
iV-dependent manner. 



C. Lack of Mean Field Effects in the 
Single-Particle Wavefunctions 

Diagonalization of the Hamiltonian (JHJ within the ba- 
sis of restricted Fock states ©, furnishes a basis repre- 
sentation of the ground state and N excited states due 
to exchanges of atoms between the two single-particle 
states \x\) and |xa)- As no equations have been derived 
for these states, their functional form is not specified by 
this approach and approximate models have been chosen 
that either define the xi and X2 as solutions of the single- 
particle Schrodinger equation or parameterize the Hamil- 
tonian matrix elements hki and Vkimn directly. Neither 
method takes into account the effects of the condensate 
mean field on the shape of the single-particle wavefunc- 
tions, which becomes more important as the interaction 
strength between the constituent atoms becomes larger. 
We will demonstrate, in Section IV, how to build atomic 
interactions into the underlying single-particle states that 



enter the matrix elements in the many-body Hamiltonian 
©. But first, it is important to discuss mean field theory 
by itself for identical bosons. 

III. REVIEW OF MEAN FIELD THEORY FOR 
IDENTICAL BOSONS 

Single-particle wavefunctions Xk were first introduced 
in the Fock space approach of the previous section. 
Within that model, no equations were developed to de- 
termine these functions. In this section, we derive a set of 
equations by minimizing the energy associated with (jSJ) 
to dictate the functional form of the two Xk ■ The equa- 
tions that arise are the bosonic Hartree-Fock equations. 

A. Restriction to One Single-Particle State 

The mean field theory in which all N atoms occupy 
the same single-particle state \x) = ^|vac) is known 
as the Gross-Pitaevskii (GP) theory. The many-body 
wavefunction ^f N restricted to the Fock state \N) = 
(w) N \vac) / >/N\ is approximated by the product 

vf H (l,...,A0 = x(l)x(2)---x(A0 (15) 

of single-particle wavefunctions x- This particular type of 
product is also called a Hartree product since it involves 
no symmetrization whatsoever. Variation of the expec- 
tation value of the many-body Hamiltonian Q with 
if! = xb m \^ H ) = \N) with respect to % and subject 
to the constraint that x is normalized to unity leads to 
the GP equation 

{{-h 2 /2m)V 2 + V ext + g(N - l)|x| 2 }x = MX (16) 

provided that the contact interaction approximation has 
been made. The chemical potential /i enters as 
a Lagrange multiplier which enforces the normalization 
(x|x) = 1- While this approximation provides an appro- 
priate mean field description of the simple BEC, it is 
not flexible enough to characterize the break up of a sin- 
gle condensate into multiple fragments, where potentially 
several single-particle wavefunctions are macroscopically 
occupied. 

B. Restriction to Two Single-Particle States 

A mean field theory, which generalizes the GP (or 
Hartree) ansatz by adding a second single-particle state, 
is the bosonic Hartree-Fock (BHF) theory. The BHF 
ansatz rests on approximating the many-body wavefunc- 
tion for N bosons as a symmetric product of the two 
single-particle wavefunctions xi and Xi 0- IIH • That is 

* BHF (l,...,iV) 

= 5{ X i(l) ■ • • Xi(N[) X2 (N[ + 1) •• -X2(N[ + N£}, 

(17) 
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where S — (vQvT) -1 J2p P is the symmetrization opera- 
tor and P is an operator that permutes the atomic coor- 
dinates. We place primes on the occupation numbers N[ 
and N 2 for reasons that will become evident in Section 
IV. This BHF wavcfunction is also called a single per- 
manental wavefunction in contrast to the single determi- 
nantal wavcfunction for fermions built from an antisym- 
metric product of single-particle wavefunctions. If we 
had omitted the symmetrization operator S all together 
in l|17fl . then we would have a two-single-particle state 
Hartree or GP theory 13] ■ Two-state Hartree theory 
provides an alternative mean field theory that neglects 
the quantum-mechanical exchange interaction associated 
with identical particles. We note that the ansatz Ijl7(l has 
been explored in a different context in Q. 

The coupled BHF equations may be determined by 
variation of the expectation value of the functional 

K\xi,X2} = H- VkiN' k (A kl - 6 kl ) (18) 

fcZ=l,2 

in the BHF state 
|* BHF ) = \N[,N$ = (St)^(^)^|vac)/V^M! (19) 

with respect to the two wavefunctions xi and X2, where 
H is the many-body Hamiltonian © restricted to the 
model space of two single-particle states, N = N[ + N 2 
is the total particle number operator, and A« = (xk\xi) 
are the matrix elements of the wavefunction overlap A. 
The second term on the right hand side of this equation 
adds Lagrange multipliers fiki whose purpose is to con- 
strain the single-particle wavefunctions to be orthonor- 
mal. This leads to the coupled two-single-particle state 
BHF equations 

hxi + (N[ - l)r lX i + N' 2 [J 2 + K 2 ]xi = M11X1 + M12X2 

hx2 + {N' 2 - l)r 2 X2 + N[[Ji + &i]X2 = M21X1 + M22X2, 

(20) 

where /i(x) = (—h 2 /2m)V 2 + V^t^x) is the sum of ki- 
netic energy and external trapping potential, T k accounts 
for the interaction of one atom in the fcth single-particle 
state with the mean field of N' k — 1 other atoms in the 
same state, Ji is the direct interaction between a sin- 
gle atom in the fcth single-particle state and the mean 
field of N! atoms in the Zth (I ^ k) single-particle state, 
and ICi is the exchange interaction between states k and I 
which arises due to the symmetrization of the BHF wave- 
function. These equations have already been derived by 
others 0, 0, |3i| and have been further extended to treat 
identical bosons in arbitrarily many single-particle states 
0- Similar equations have been for derived for distin- 
guishable multicomponent (spinor) BECs [28l I39L |40L Elf . 
however, it is important to note that they are not the 
same as the BHF equations <|2U|) . which describe identi- 
cal bosons. 

The diagonal Lagrange multipliers [ikk in (|18(l ensure 
the proper normalization of the single-particle wavefunc- 
tions while the off-diagonal fiki enforce their orthogonal- 
ity. If the external potential V C xt is symmetric, then the 



off-diagonal Lagrange multipliers are not necessary as the 
wavefunctions are automatically spatially orthogonal by 
symmetry. The fiki cannot in general be removed by uni- 
tary transformation as in the fermionic case 0, |(| . Con- 
sequently, arbitrary linear combinations of xi and X2 are 
not solutions of the BHF equations (|20|> . Koopmans' the- 
orem IUI23, E3 is satisfied for the diagonal multipliers. 
Therefore //u = E BBF [N[,N^} - E BBF [N[ - 1, A^] and 
A*22 = E BiiF [N[, N!,} - E BliF [N[, Nl, - 1] inherit the roles 
of chemical potentials 0, 0|, where the BHF energy is 
given by 

E BBF = Khkk + (1/2) N k( N 'k ~ Wkkkk 

k=l,2 k=l,2 

+ (1/2) KNliVkiki+VkUk], 

fc^Z=l,2 

(21) 

and where hu and Vki mn are matrix elements of the 
single-particle Hamiltonian h(x) and two-body interac- 
tion potential V(x, x') in JJJ). The direct and exchange 
integrals in (|20|l are defined as 

r fe (x) Xfc (x) = / v [ Xfc (x')F(x,x')xfe(x')]Xfc(x)d 3 x' 
JKx)Xfc(x) = / ^/ [x ^ (x')^(x,x')x ^ (x')]Xfe(x)d 3 a; , (22) 
K,(x)xfc(x) = / v [ x Kx')^(x,x')x fc (x')]Xi(x)rf 3 x' 

for k 7^ I and k,l — 1,2. Note that the potential Tk is 
a direct interaction that arises only for bosons. There is 
no analogous term for fermions due to Pauli exclusion. 

The BHF theory reduces to the GP theory in the ex- 
treme limit of only one occupied single-particle state, i.e., 
when N[ = N and N 2 = or vice versa as demonstrated 
in Figure 13 However, unlike the GP equation l|16(l . the 
BHF equations H2U|) offer additional freedom in that they 
accommodate the possibility of various numbers of atoms 
in each of two single-particle states. By choosing N[ 
atoms to be in the first state and N 2 atoms to be in the 
second, a set of BHF equations corresponding to that 
particular arrangement is obtained. Each arrangement is 
associated with a symmetrized BHF single-permancntal 
wavefunction like l|17f) . Where BHF equations specify 
the underlying single-particle wavefunctions xi and X2-, 
we call the Fock state |* BHF ) = \N[,N^) a single con- 
figuration. These Fock states are distinguished from 
those in Section II, which were written without primes 
as \Ni,N 2 ), because the single-particle states that make 
up each single configuration are determined by solving 
BHF equations. It is then possible to vary the number of 
particles in each of the single-particle states to find the 
lowest energy single configuration for a particular trap 
geometry. This flexibility gives rise to two very different 
physical meanings for xi and X2- In the first, N[ ~ N 
and there are a relatively small number of atoms in |x2)- 
In this regime, |xi) is a condensate state and IX2) is a 
single-particle excited state. The energetic cost of mak- 
ing excitations from |xi) to IX2) is macroscopically large 
due to the effect of bosonic amplification In the 



7 




N'JN 



FIG. 3: (Color online) The BHF energy reduces to that of 
the GP ground state when N[/N = 1 and to that of the 
first GP excited state when N[/N — 0. This corresponds to 
all particles in the symmetric GP ground state and the first 
antisymmetric GP excited state respectively. A dimensionless 
interaction strength of ojqid = 40 [see lUTl l for fOO atoms m 
a single well potential was used in this calculation of E BHF . 



second, fragmented, regime, N[ and are both on the 
order of N. In this case, both single-particle states are 
condensate states. Single-atom excitations between 
and IX2) also cost a macroscopic amount of energy, but 
now there is also a macroscopic gain of energy as the 
atom is added to a second macroscopically occupied state 
|l8j . This distinction, which is quite important, will be 
elaborated on in Section V. 



C. Numerical Integration Method 

We numerically solve the coupled BHF equations us- 
ing a fast Fourier transform based pseudospectral grid 
method |43|, |44J in quasi-one dimension |45|, |46j. Note 
that quasi-one dimension does not mean one dimension 
but rather that the variation of the BEC density is neg- 
ligible in the two transverse dimensions and a separa- 
tion of variables is permissible. Using this approach, 
the equations are expanded onto a discrete Fourier sine 
basis with 2 s fixed grid points j^, which satisfies the 
proper boundary conditions, and the expansion coeffi- 
cients are variationally optimized. Rather than solving 
the time-independent equations (|2U[1 self-consistently, we 
work with their time-dependent version, where /in and 
/i22 are replaced by ih(d/dt) j2^. We then solve the 
time-dependent BHF equations by the method of steep- 
est descents in imaginary time |27| . That is, we employ 
the Wick time rotation t — > t = it, which takes ih(d/dt) 
to —h(d/dT), and integrate the pair of coupled nonlinear 
diffusion equations 



as an initial value problem in r, where J~k is the boson 
Fock operator for the kth single-particle state (k = 1,2), 

^■ fc (x) = ^(x) + (^-I)r fe (x) + 7V/[J / (x) + /C / (x)]. (24) 

The resulting equations are well-defined once the two- 
body interaction potential V is identified. As was done 
in Section II, we make use of the contact potential 
V(x, x') = (4:Trh 2 a/m)5(x — x'). Following the argument 
of Esry et. al. |2q , the renormalized S- wave scattering 
length a is taken from multichannel T-matrix calcula- 
tions using symmetrized two-body wavefunctions. Thus, 
the contact interaction approximation effects the replace- 
ment of the symmetric combination of two-body matrix 
elements in the BHF energy Ij21|l with a single contact 
potential and not each matrix element separately. In 
symbols that is 

Vhm + V kUk = (kl\V\kl + Ik) -» (4nh 2 a/m,)8( X - x'). 

25 



-H(d/dT)xi = T1X1 - M12X2 
-h{d/dr)x2 = F2X2 ~ M21X1 



(23) 



This identification differs by a factor of two from 
Q; see [Sij, [EH anc l ^he discussion following (|39|l . The 
resulting BHF equations (H^ are [HHH 

{h + g(N[ - l)\xi\ 2 + gN^\ X2 \ 2 }xi = M11X1 + 

{h + g(Nl, - l)\x 2 \' 2 + gN[\ X i\ 2 }X2 = A*2iXi + ^22X2- 

(26) 

We choose as an initial condition for the single-particle 
wavefunction xi the square root of the Thomas-Fermi 
density p TF (x) = \p, — V cx t(x)]/g{N — 1) at a particu- 
lar barrier height, chemical potential, and value of cou- 
pling constant. For a symmetric trapping potential, the 
second wavefunction xi is taken to be antisymmetric 
to xi and the off-diagonal Lagrange multipliers are un- 
necessary. Otherwise, xi need only be orthogonal to 
Xi but constraints are needed to maintain orthogonal- 
ity. Both wavefunctions are initially normalized so that 
llxi(0)l! = ||X2(0)|| = I and are real-valued. Thus, the 
overlap matrix elements A^; = S^i initially. We then em- 
ploy the standard relaxation approach, i.e., we subtract 
a guess fL of the ground state chemical potential from the 
Fock operator and allow the system to time evolve. The 
wavefunctions relax according to 

X fc (r) « 5>xp(-[ M - fi kk ]r/h)xt(0K (27) 

for k — 1,2, where the expansion coefficients c M are pro- 
jections of the evolving state |Xfe( r )) onto the station- 
ary basis \Xk( T ))- Eventually, all excited states decay 
away after repetition of this procedure together with in- 
termittent renormalization. If the wavefunctions are ini- 
tially orthogonal and share the symmetry of the trapping 
potential, then the wavefunctions that remain are the 
symmetric and antisymmetric solutions which minimize 
_B BHF for each configuration. 

Whenever the external potential V e xt is asymmetric, 
off-diagonal Lagrange multipliers must be introduced to 
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keep the single-particle wavefunctions orthogonal to each 
other throughout the evolution. The proper choice for 
the /iki (k 7^ I) will ensure that {d/ dr)Aki(nki) = for 
all time t. The appropriate fiki (k ^ I) are found by 
multiplying the first equation in 123(1 by X2 and the sec- 
ond equation in (|23jl by xi, adding the two equations to- 
gether, and then performing a volume integral. One then 
arrives at the following expression for the time derivative 
of the off-diagonal overlap 



-H{d/dT)A 12 = (X2\Fl - MXllXl) - Ml2(X2|X2) 

+ (Xl\^2 ~ M22IX2) - M2i(Xl|Xl), 



(28) 



where the Jlkk are intermediate guesses of the ground 
state chemical potentials associated with xi and X2 ■ The 
time derivative of A21 is the same, as A T = A. With 
(|2"5|) and the relation N[/j, 12 = N 2 fi 2 i @, values for \i\ 2 
and H21 can be chosen so that the right hand side of l|28|) 
equals zero. Those values are 



A* 12 



M21 



{X2\F\ ~ /fnjxi) + (xil-^2 " M22IX2) 
(JV(/J^)(xi|xi> + <Xa|X2) 

{X2\F\ - gnjxi) + (xil-^2 ~ M22IX2) 
(XilXi) + {NUK){X2\X2) 



(29) 



With this choice, the time derivatives of the off-diagonal 
matrix elements of A are zero. Therefore, if the single- 
particle wavefunctions are initially orthogonal, then they 
will stay orthogonal for all time r regardless of trap sym- 
metry. For each choice of BHF configuration \N[,N£), 
the relaxation algorithm will then find the associated 
two orthogonal wavefunctions which minimize the energy 
l|21jl. This is the essence of our integration scheme. 



D. Hartree-Fock and Schrodinger Wavefunctions 

Figure 0] displays the BHF single-particle wavefunc- 
tions xi and X2, which are obtained by solving l(26jl in 
quasi-one dimension at four different barrier heights 0, 3, 
9, 13 hut. In each panel, Xi is associated with the BHF 
configuration \N, 0) and x.2 is associated with the BHF 
configuration |0, TV). In these extreme configurations, 
the single-particle BHF wavefunctions reduce to the GP 
ground and first excited wavefunctions. The configura- 
tion \N, 0) minimizes E BUF at each barrier height. Figure 
displays the solutions of the single-particle Schrodinger 
equation hxk — £kXk at the same four barrier heights 
0, 3, 9, 13 hw. Since the Schrodinger equation includes 
no atom-atom interaction, the Schrodinger solutions have 
no dependence on the number of particles in each single- 
particle state. For 100 atoms with a dimensionless inter- 
action strength of ckqid = 40 [see (|47f) ]. it is seen that 
the BHF wavefunctions follow the Thomas- Fermi result 
quite well whereas the Schrodinger wavefunctions bear 
little resemblance to either BHF or Thomas-Fermi solu- 
tions. 



E. Lack of Correlation Effects 

The solutions xi and %2 of the BHF equations (|20|l in- 
clude the effects of the condensate mean field on their 
shape. An appropriate symmetrized product of these 
wavefunctions results in a BHF wavefunction like that in 
p7[l. Unlike in the single determinantal case where the 
wavefunction is defined only up to unitary transforma- 
tions of its constituent single-particle wavefunctions ^3 > 
the BHF single permanental wavefunction seemingly cor- 
responds uniquely to a single configuration. That is, the 
single permanent 

5{xi(l) ' ' -Xi(N[)x2(N[ + 1) • • • X2(N[ + N^)} (30) 

is in one-to-one correspondence with the single configu- 
ration \N[, N2), while another BHF permanent 

S{xi(X) ■ ■ ■ Xi{NDxi{N'{ + 1) ■ • • X2(N[' + N 2 )} (31) 

is in one-to-one correspondence with the single configura- 
tion \N", N 2 ), and so on. All such single configurations 
can be collected into the set 

{\N, 0), \N - 1, 1), |JV - 2, 2), . . . , |0, N)}, (32) 

which, like ©, spans an (N + l)-dimensional Fock space 
but, in addition, has states that are now specified by BHF 
equations. 

Considering that the eigenstates of the many-body 
Hamiltonian © are not, in general, a single configura- 
tion but are rather composed of a linear combination of 
N + 1 such configurations, it is evident that a single con- 
figurational description is quite limiting. In particular, it 
lacks all effects of correlation that arise between atoms 
in different configurations. To this end, we formulate a 
new approach in Section IV that combines the full diago- 
nalization of the many-body Hamiltonian in a restricted 
basis with mean field theory for the underlying single 
particle states. 



IV. MULTICONFIGURATIONAL BOSONIC 
HARTREE-FOCK THEORY 

The many-boson theory restricted to a Fock basis con- 
sisting of two single-particle states and the two-single- 
particle state BHF mean field theory provide comple- 
mentary descriptions of the BEC and its fragmentation 
into two condensates. For arbitrary interaction strength, 
the BHF approach is well justified whenever the state of 
the BEC can be described by a single configuration, but 
breaks down whenever a multiconfigurational description 
is appropriate. Alternatively, the finite basis represen- 
tation of the many-body Schrodinger equation accounts 
for atomic correlation between each Fock state and in- 
cludes the effects of the condensate mean field directly in 
the Hamiltonian. However, it does not provide equations 
that specify the underlying single-particle wavefunctions. 
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FIG. 4: (Color online) BHF single-particle wavefunctions xi an d X2 versus coordinate z at four different barrier heights 0, 3, 
9, 13 huj. A dimensionless interaction strength of oqid = 40 [see 147H 1 was used. In each panel, xi corresponds to the BHF 
configuration |iV, 0), while \2 corresponds to the BHF configuration |0, iV). The BHF energies associated with |iV, 0) and |0, N) 
at zero barrier height appear in Figure Both wavefunctions have been set at the chemical potentials /in and (122 associated 
with each configuration. At each barrier height, we plot the corresponding Thomas-Fermi wavefunction as a solid black curve. 



The use of parameters or even single-particle Schrodinger 
wavefunctions may not capture certain properties of the 
BEC that are associated with strongly interacting atomic 
gases, where mean field effects on the shape of the single- 
particle wavefunctions are important. 

To this end, we variationally combine the BHF mean 
field theory of Section III with the restricted Fock state 
representation of the many-body theory of Section II, 
allowing for both: 

(1.) The effects of the condensate mean field on the 
shape of the single-particle wavefunctions. 

(2.) The ability to describe states that are made up of 
multiple configurations. 

The former is necessary in the strongly interacting 
regime, while the latter is needed to describe conden- 
sate fragmentation within our approach. We refer to the 
union of these disjoint theories as the multiconfigurational 



bosonic Hartree-Fock theory or MCBHF. This theory is 
rich enough to characterize the atomic structure of the 
simple BEC and its fragmentation into two condensates 
at zero temperature as its trapping potential is deformed 
from a single well to a double well with large barrier 
height. 



A. General Theory 

The basic idea behind our MCBHF approach is to diag- 
onalize a representation of the many-body Hamiltonian 
© in a set of basis functions of the form 

\Ni,N 2 ',{N{,N£), (33) 

where the total number of atoms N = Ni+N 2 = N[+N%, 
and N±,N{ — 0, . . . , N. These basis states are a com- 
bination of the Fock states of Section II and the BHF 
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FIG. 5: (Color online) Schrodinger single-particle wavefunctions xi an d X2 versus coordinate z at four different barrier heights 
0, 3, 9, 13 hit). These solutions have no dependence on the number of atoms in each single-particle state. Both wavefunctions 
have been set at the same BHF chemical potentials used in Figure 0] At each barrier height, we plot the Thomas-Fermi 
wavefunctions of Figure[I]as a solid black curve. Little resemblance is seen between the Schrodinger and Thomas-Fermi results. 



configurations of Section III. It will become evident that 
iVi and N2 count atoms in left- and right- localized Fock 
states, while N[ and count atoms in symmetric and 
antisymmetric BHF states. We call the kets gen- 
eralized configuration states or GCSs because for every 
underlying BHF configuration \N[, N^), which has been 



indicated in by {N[,NZ,}, there are an additional 
N + 1 Fock states |iVi, JV 2 ) that can be built from this 
BHF reference configuration. For instance, the subset of 
GCSs stemming from the particular BHF configuration 
\Ni,N$ is 



{\0,N;{N[,N£), \l,N-l;{N[,Nfi), \2, N - 2; {N{, N^}, |iV,0;{iV{,^})}. 



(34) 
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However, there are (TV 
organized into the set 



l)-many underlying BHF configurations in total. The collection of all GCSs (|33[1 can be 



|0,iV;{0,iV}), 
\Q,N;{1,N- 



\0,N;{N{,N£), 



1,7V- 
1,7V- 
1,7V- 


1;{o,tv}), 

1;{1,JV-1», 
l;{2,7V-2}), 


|2,7V- 
|2,7V- 
\2,N- 


2;{0,iV}), 
2;{1,7V-1}), .. 
2;{2,7V-2}), .. 


., |7V,0;{0,7V}), 

|7V,0;{1,7V-1}>, 
., |7V,0;{2,7V-2}), 


1,7V- 


l;{N[,N^}), 


\2,N- 


2;{N{,N>}), .. 


., |7V,0;{7V(,7Va>, 


1,7V- 


i;W0}>, 


12,7V- 


2;{7V,0}), 


., |7V,0;{7V,0}) 



(35) 



I |0,7V;{7V,0}), 

containing (TV + 1) x (TV + l)-many elements of which (|34|l makes up just one row. These elements span a Fock space 
that contains the (TV + l)-dimensional Fock spaces of Sections II and III in certain limits. For example, the set of all 
GCSs (|35J) reduces to the set (01 associated with the Schrodinger model of Section II in the limit of vanishing atom- 
atom interaction within the BHF equations. This corresponds to the subset of all TV + 1 columns belonging to any 
single row of (|35|l . since all rows within each column are equivalent in the noninteracting limit of BHF. Alternatively, 
the subset of GCSs of the form \N 1} TV 2 ; {N[ = TVi,TV^ = TV 2 }) is identical to the set of BHF configurations ^ of 
Section III. The diagonal entries of l|35|l realize such a subset. 

For each BHF configuration |7V{,7V2), which is determined by solving the BHF equations it2U|) . the MCBHF state 
vector can be expanded onto the GCS basis (|33() according to 



N 
ATi=0 



(36) 



where the expansion coefficient C V N [TV{ , is the probability amplitude for the vth excited state of the system to be 
in the GCS |TVi, TV 2 ; {TV{, TV^}). With the MCBHF states, the many-body Schrodinger equation 



H\* N ; {N[, N£)„ = E"[N[ , N']\* N ; {N[, N^}„ 



(37) 



is represented in the GCS basis (|33|l . In a symmetric trapping potential V^ x t, the Hamiltonian matrix elements 
hki = (k\h\l) and Vkimn — (kl\V\mn) are constructed out of the left- and right-localized states (l/v2)[|xi) ± \x2)], 
which are obtained from the symmetric and antisymmetric BHF solutions \xi) and |% 2 ) by unitary transformation. 
Such a transformation can always be performed within h^i and Vkimn since they provide only matrix elements for the 
many-body Hamiltonian 10. Consequently, the Fock space occupation numbers N\ and 7V 2 , which are the first two 
entries of the GCS l|33|l . refer to the number of atoms in left- and right-localized single-particle states, while the BHF 
occupations numbers 7V{ and TVj, which appear in brackets within (|33|l . refer to the number of atoms in symmetric 
and antisymmetric single-particle states. Full diagonalization in this basis gives rise to a set of energy eigenvalues 
E v [TV{ , TVj] and associated eigenvectors C v Ni [TV{ , TV^] that depend upon the particular BHF reference configuration 
|TV{, TV2). Solving this eigenvalue problem for each of the 7V + 1 configurations results in the set of all MCBHF energies 



{E v [0,N\, E"[1,N-1], E v %N-2], ^[TV{,TV 2 ], 



r 



E v [N,0]}, 



(38) 



, TV. The cardinality of this set is (TV 



where v = 0, 
1) x (7V + 1). 

We then appeal to the variation al p rinciple |48| and 
Hylleraas-Undheim theorem |49ll5Ctl5lll52| . which is pre- 
sented in the Appendix without proof. The variational 
principle asserts that within the model space of MCBHF, 
the ground state energy E° stemming from any BHF con- 
figuration is an upper bound to the exact ground state en- 
ergy and may be systematically reduced by adding more 
GCSs. Furthermore, by the Hylleraas-Undheim theorem, 
we know that ^th MCBHF excited state energy E v stem- 
ming from any BHF configuration is an upper bound to 



the exact i/th excited state energy, and may also be sys- 
tematically reduced by adding more GCSs to the basis. 
Therefore, it is true in general, that the optimal BHF 
configuration at each energy level is the one that mini- 
mizes the energy at that particular level. By choosing the 
BHF configuration |TV{,TV2) that minimizes E°, and the 
BHF configuration ITV", TV2'} that minimizes E 1 , and so 
on, we generate a set of TV + 1 optimal MCBHF energies 



{E°[N' 1: 



N' 2 ] 



E 



[N", TV 2 ] 



E»[N[",N?]}, (39) 
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where each energy E v may stem from a different BHF 
reference configuration. In fact, for each configuration, 
the coupling constant g in the BHF equations (|26|1 could 
also be treated as a variational parameter [61], but in this 
case g may not be the same as the physical coupling con- 
stant. Each of these optimal energies has an associated 
optimal eigenvector that is based off of the same BHF 
configuration that optimizes the energy. For example, 
the set of eigenvectors corresponding to the set (|39|) is 

{C° Nl [Ni,NA, C l Nl [N{',N^}, C^[Ni",N^'}}, (40) 

where iVi = 0, . . . , N and where each of the individ- 
ual BHF reference configurations {N{,N£}, {N{',N%}, 
and {N[",N%'} are the same in both and gQJ. In 
this way, the MCBHF theory simultaneously optimizes 
both the underlying basis functions, i.e., the BHF ref- 
erence configurations, and the expansion coefficients of 
the GCSs. Therefore, each state vector \^ N ; {N[, N^})^ 
is potentially made up of multiple GCSs between which 
arise correlations, and where the effects of the conden- 
sate's mean field are built into the Hamiltonian as well 
as into the shape of the underlying single-particle wave- 
functions. 

We call the MCBHF state vector that has the lowest 
MCBHF energy at the vth energy level the optimal i/th 
MCBHF state vector and denote it by 

|**5 «,<}>„. (41) 

The bullets symbolize that the energy of the i/th state 
is minimized in the BHF configuration |iV{ , iV^}. The 
corresponding set of constituent GCSs forms an opti- 
mal (N + l)-dimensional subset of (|35(l that spans a 
Fock space supporting lower energies than either of those 
in Sections II and III. We now provide an algorithmic 
method by which to realize this subset. 



B. Implementation Algorithm 

Implementation of the MCBHF theory begins by solv- 
ing the BHF equations l|26|l for a particular trap geom- 
etry, say a single well potential, and a particular BHF 
configuration, say the single configuration j6^ 

\N{ =0,N^=N) = |* BHF ) = (b\)°(bl) N \vac)/V0\Ni., 

where there are no atoms in the single-particle state 
and N atoms in the single-particle state \x2}- We choose 
V cx t to be symmetric so that the solutions xi an d Xi 
of the BHF equations are symmetric and antisymmet- 
ric. With these single-particle wavefunctions we first cal- 
culate the BHF energy (|21|l . We then unitarily trans- 
form xi a- n d Xi to the left- and right-localized functions 
(l/\/2)[xi ± X2L and use them to build the matrix ele- 
ments hki and Vkimn in the many-body Hamiltonian JfJJl . 



With the MCBHF state vector 

N 

\* N ;{0,N}) V = ]T C& l [0,JV]|JV 1 ,JV-iVi;{0,JV}>, 

7Vi=0 

full diagonalization of the Hamiltonian yields a basis rep- 
resentation of the ground and excited state energies and 
associated eigenvectors stemming from the single BHF 
configuration \N[ = 0, Nl, = N). 

Next, we solve the BHF equations <|26[) associated 
with the single configuration \N[ = 1, N% = N — 1) = 
(&J) 1 (S|) JV_1 |vac)/ A /l!(iV — 1)!, in which a single atom 
has been moved from the second single-particle state to 
the first. With the BHF solution for this configuration we 
again calculate i? BHF , make a unitary transformation to 
left- and right-localized states, and then build the matrix 
elements hu and Vkimn- With the MCBHF state vector 

1*^(1,^- 

JV 

= C^N-lWN^N-N^ihN-l}) 

N 1= 

(44) 

the full diagonalization of the associated Hamiltonian 
yields the MCBHF ground and excited state energies and 
eigenvectors stemming from \N[ = 1, iV^ = N — 1). 

We then repeat this procedure for each BHF con- 
figuration until we reach the last BHF configuration 
\N,0) = (b\) N (bl)°\vac)/y/NW. in which all N atoms 
are in the symmetric state \xi) and no atoms are in 
the antisymmetric state \x2}- Mapping out the complete 
set of MCBHF energies versus all possible BHF refer- 
ence configurations at this trap geometry results in a 
diagram like that displayed in Figure HO Notice the 
detailed energy level mergings and splittings that oc- 
cur in this figure. By appealing to the variational and 
Hylleraas-Undheim theorems, we find the variationally 
optimal MCBHF ground state and excited state energies 
and denote their location by a single bullet. The location 
of these bullets also indicates which BHF configuration 
provides the optimal MCBHF ground and excited state 
energy. In this way, we are variationally optimizing both 
the shape of the single-particle wavefunctions and the 
expansion coefficients. That is, a nonlinear optimization 
determines the shape of the BHF single-particle wave- 
functions for each BHF configuration while a second lin- 
ear Hylleraas-Undheim optimization determines the ex- 
pansion coefficients of the optimal MCBHF state vector 
{N{',N£}) V . We then repeat this algorithm as the 
trapping potential is deformed from a single to a double 
well geometry. 

C. Example: MCBHF Excited State 

Consider, for example, the optimal 51st excited state 
in Figure© It is built off of the \N{ = 56, = 44) BHF 
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FIG. 6: (Color online) MCBHF ground and excited state ener- 
gies and BHF ground state energy for N = 100 atoms in a sin- 
gle well potential (or zero barrier limit of a double well) with 
a dimensionless coupling of qqid = 40 [see 1471 1 qqid = 40 
[see 1471 1 versus BHF reference configuration. The solid black 
curve represents the BHF ground state energy. Bullets are 
placed at the minima of BHF ground and MCBHF ground 
and excited energies. The line style is alternated to aid in 
visualization. This figure is a zoomed in and more detailed 
version of the first panel in Figure |5] 




FIG. 7: (Color online) MCBHF expansion coefficients Cffc 
of the 51st excited state \^ N ; {56*, 44*})si versus the occu- 
pation number in each CCS. The label {56,44} in the state 
vector signifies that the minimal 51st MCBHF excited state 
energy has been attained in the BHF configuration \N[ = 
56, N2 = 44) . This probability amplitude corresponds to 100 
atoms in a single well trapping potential with oqid = 40 [see 



reference configuration, which can be read off of Figure 
[SJ and has expansion coefficients Cjy [56, 44] as shown 
in Figured This particular optimal MCBHF excited 
state is not described by the single BHF configuration 
\N[ = 56, N£ = 44) displayed in Figure |S] but rather by 
almost all possible GCSs within this model space. It is 
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FIG. 8: (Color online) BHF single-particle wavefunctions \i 
and X2 versus coordinate z associated with the BHF con- 
figuration \N{ = 56, N2 = 44), which underlies the optimal 
51st MCBHF excited state at a barrier height of Hui. A di- 
mensionless coupling of «qid = 40 [see 1471 1 for 100 atoms 
was used. Both wavefunctions have been set at the chemi- 
cal potentials fin and ^22 associated with this configuration. 
The corresponding Thomas-Fermi wavefunction is plotted as 
a solid black curve. 



<H3] 



roughly of the form 



* JV ;{56*,44 , }) 51 sa <7f 1 |10,90;{56,44}) + Cf 1 |ll,89;{56,44}) + ---^ 

(45) 



where Cjy = C]y i [56,44] and where the label {56,44} many-body Hamiltonian © for the optimal 51st excited 
in each ket signifies that the matrix elements of the 
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state have been built out of the single BHF configura- 
tion \N[ = 56, = 44), where N = 100 = iV x + N 2 . 
While this BHF configuration minimizes the energy for 
the 51st excited state, another BHF configuration may 
be optimal for the 52nd, and so on. Furthermore, each 
optimal state may also have a different distribution of ex- 
pansion coefficients Cjy ■ In other words, every optimized 
MCBHF state is potentially derived from a different BHF 
reference configuration and is potentially made up of a 
unique superposition of GCSs. 



V. RESULTS AND DISCUSSION 

We employ an external trapping potential V ex t that 
is harmonic with a Gaussian function centered at the 
oscillator's minimum. In quasi-one dimension |45l l46| . it 
has the dimensionless form 

V^xtO) - (l/2)[z-4/2] 2 + Aexp(-[z-4/2] 2 /2), (46) 

where t z is the spatial grid length and A is the ampli- 
tude of the Gaussian having unit width in oscillator units. 
While we do not intend to precisely model a particular 
experiment in this paper, the functional form of V ex ± |53j ] 
has been chosen because of its resemblance to the double 
well interference experiments performed at MIT 0,0)13 • 
In practice, A — 1 + eB and B is varied from zero up to 
150 in increments of e = 0.1 fuo. As the amplitude is 
increased, the constant 1 + In A is subtracted off of the 
potential so that its minimum always lies at hu. All 
energies are, therefore, reported as relative to the zero 
of the trapping potential. In this way, V ey ± is deformed 
from an initially single well with zero barrier height to a 
double well with a barrier height of eB = 15 Hui. 

In the GP equation (|16fl . varying g is the same as vary- 
ing N. However, g and N enter the many-body Hamilto- 
nian © in different ways so that varying g is not the 
same as varying N in MCBHF theory. Although ac- 
tual experiments involve N w 10 5 condensate atoms, we 
choose to perform all calculations with a smaller number. 
Nevertheless, it is possible to illustrate the importance 
of correlation and mean field effects by adjusting g so 
that the product gN has the correct value. In quasi-one- 
dimension, the dimensionless coupling constant becomes 

agio = inaP z N/Ll, (47) 

where (3 Z = \Jh/muj z is the z-oscillator length and L± 
is the transverse length associated with the quasi-one- 
dimensional approximation. Choosing trap frequency 
uj z = 2tt x 30 Hz and S'-wave scattering length a = 4.9 
nm appropriate for the recent 23 Na double well inter- 
ference experiments 0, , the dimensionless quasi-one- 
dimensional coupling constant has an approximate value 
of 50, where a condensate density of 10 13 atoms/(cm) 3 
has been assumed. In light of this value, we perform 



MCBHF calculations at oiqid = 40 where Thomas- 
Fermi characteristics are already evident. The quasi- 
one-dimcnsional BHF wavefunctions displayed in Figure 
^correspond to this particular value of coupling strength. 
Throughout the paper, we have chosen the frequency 
scale lo to be lo z . 



A. MCBHF energy versus BHF configuration 

Allowing the intcrwell barrier in (|46|l to grow, we re- 
peat the implementation algorithm of Section IV B for 
each and every barrier height. A panorama of MCBHF 
energies versus BHF configuration is presented in Fig- 
ure corresponding to four different barrier heights 0, 3, 
9, and 13 Hlj with Q!qid = 40 in both the many-body 
Hamiltonian © and in the BHF equations (|26|l . The 
first panel is an enlarged version of Figure where only 
every fifth line is plotted to aid in visualization. Bullets 
have again been placed at each of the variationally opti- 
mal energies. The structure of the MCBHF ground and 
excited state energies versus BHF reference configuration 
is quite detailed and it is important to explain some of 
its features. For example, it can be seen in FigureEJor in 
the upper two panels of Figure^that there are two ener- 
getic pathways that can support optimal solutions. One 
pathway provides the globally optimal solutions while the 
second pathway provides only locally optimal solutions. 
One might ask why this structure is present and why it 
is not symmetric about the \N/2, N/2) BHF configura- 
tion. The answer lies in the role of the two single-particle 
states \xi) and \x2), which are the symmetric and anti- 
symmetric BHF solutions, as various numbers of atoms 
are placed in each state. 

Consider the BHF reference configuration |0, N) in 
which there are no atoms in the symmetric state and 
N atoms in the antisymmetric state. Excited MCBHF 
states |\t ; {0,N}) U are made up of single atom excita- 
tions out of the highest energy BHF reference configu- 
ration describable in the basis. Furthermore, these exci- 
tations move atoms to states that have zero occupation. 
Due to the effect of bosonic amplification, these excita- 
tions of an excited state cost a macroscopic amount of 
energy [TT| . Alternatively, consider the MCBHF state 
^^{A^O})^ which is built out of the BHF reference 
configuration |./V, 0), in which all atoms are in the sym- 
metric ground state. Single atom excitations out of this 
state also cost a macroscopic amount of energy, due 
to bosonic amplification, as atoms are moved from a 
macroscopically filled state to empty states. However, 
these excitation energies lie lower than the previous en- 
ergies because they correspond to excitations out of the 
ground state rather than excitations from an excited 
state. This explains why the energies on the left side of 
each panel, i.e., N[/N ~ 0, are larger than those on the 
right, i.e., N[/N ~ 1. Lastly, consider the MCBHF state 
l*^;^ w N/2,Nt, w N/2}) v , which is approximately 
built out of the BHF reference configuration \N/2, N/2), 
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FIG. 9: (Color online) MCBHF ground and excited state energies and BHF ground state energy for oqid = 40 and N = 100 
versus BHF configuration. Four different barrier heights 0, 3, 9, and 13 hu> are displayed from left to right and top to 
bottom. The solid black curve in each panel represents the BHF ground state energy at each barrier height and for each BHF 
configuration. Bullets are placed at the minima of the BHF ground state and all MCBHF energies. These are the variationally 
optimal states at each energy level. To aid in visualization, we alternate line style and only display every fifth excited state. 
The box in the first panel corresponds to Figure |S| 



in which both states are macroscopically occupied. Sin- 
gle atom excitations out of \\S? N ; {N/2,N/2}) V exchange 
atoms between two macroscopically occupied states. In 
distinction to the previous two cases, here there is a 
macroscopic energy cost to move a single atom out of 
|xi), but there is also a macroscopic energy gain as the 
atom is added to \xa)- Therefore, the MCBHF energies 
in this region represent the smallest energy excitations 
possible within our model. The location of each mini- 
mum is biased to right because it is energetically easier 
to make excitations out of the ground state than out of 
excited states. This is why the energy pathway on the 
right provides the global minimum solutions while the 
energy pathway on the left supports only local minima. 



B. Optimal MCBHF energy versus barrier height 



By collecting the set of all optimal energies at each bar- 
rier height, i.e., by collecting all the bulleted energies ap- 
pearing in Figure [5] plus those at all other barrier heights, 
an energy level correlation diagram similar to those in 
Figures n an d 121 can be made. The MCBHF correlation 
diagram associated with Figure El is displayed in Figure 
ITU1 There are qualitative similarities to those of the pa- 
rameterized and Schrodinger model calculations such as, 
for example, the pronounced ridge structure that marks 
the merging of nondegenerate energy levels to degenerate 
as the interwell barrier is raised. However, large quantita- 
tive differences are apparent by comparing the MCBHF 
correlation diagram in Figure 1101 with the correspond- 
ing Schrodinger correlation diagram displayed in Figure 
1111 The latter figures are computed by diagonalizing the 
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FIG. 10: (Color online) MCBHF energy level correlation diagram versus barrier height for qqid = 40 and A*' = 100 atoms. 
Every fifth energy level is plotted in the left panel and the line styles are alternated in both to aid in visualization. The solid 
black curve in the left panel denotes the BHF ground state energy. Energy level mergings, which lead to a pronounced ridge 
structure separating nondegenerate oscillator like states from doubly degenerate macroscopic superposition states, are focused 
on in the second panel. This panel is a zoomed and more detailed version of the boxed region in the first, where every line has 
been displayed. 
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FIG. 11: (Color online) Schrodinger energy level correlation diagram versus barrier height for oqid = 40 and N = 100 atoms. 
Every fifth energy level is plotted in the left panel and the line styles are alternated in both to aid in visualization. The second 
panel is a zoomed and more detailed version of the boxed region in the first, where every line has been displayed. 



many-body Hamiltonian © with oiqid = 40, where the 
underlying single-particle wavefunctions are taken from 
the single-particle Schrodinger equation. No optimiza- 
tion is necessary in this case since the wavefunctions have 
no dependence upon N[ and N%- To aid in visualization, 
the second panels in Figures lTUI and lTTl have been zoomed 
in precisely at the phase transition between nondegener- 
ate and degenerate energies. Notice the large difference 
in energy scales as well as the shape and placement of 
the phase transitions in Figures ITUI and ITT1 



C. Distribution of GCSs and BHF wavefunctions 

Examining the distribution of GCSs in each of the op- 
timal MCBHF states in Figure 1101 we find that below 
the ridge the distribution is harmonic oscillator like in 
form, while exotic macroscopic superpositions of GCSs 
exist above the ridge. This behavior is generic across 
all barrier heights where ridge structures are present. In 
order to exemplify these characteristics, we present op- 
timal MCBHF expansion coefficients in Figure IT21 taken 
from Figure ITUI at the particular barrier height of 3 hoj. 
Coefficients are presented that correspond to the opti- 
mal MCBHF ground state and first and second excited 
states, which are below the ridge, as well as to the opti- 
mal MCBHF 63rd, 64th, and 65th excited states, which 
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FIG. 12: (Color online) Left two panels: Expansion coefficients of the optimal MCBHF ground state and first and second 
excited states versus occupation number in each GCS as well as expansion coefficients of the optimal 63rd, 64th, and 65th 
excited states versus occupation number in each GCS. The optimal BHF reference configurations |JVi , JV2) have been provided 
for each state separately within the text. Right two panels: Expansion coefficients of the Schrodinger ground state and first 
and second excited states versus occupation number as well as expansion coefficients of the Schrodinger 63rd, 64th, and 65th 
excited states versus occupation number. These probability amplitudes correspond to N = 100 atoms at a barrier height of 3 
flu. The left column has the dimensionless coupling agio = 40 in both the many-body Hamiltonian and BHF equations, while 
the right column has oqid = 40 only in the many-body Hamiltonian. 



are above the ridge. An analogous presentation is made 
in Figure IT21 for the Schrodinger based coefficients asso- 
ciated with Figure ITTI The two left panels in Figure IT51 
correspond to the full MCBHF theory with an interac- 
tion strength of chqid = 40, while the two right panels 
correspond to the case where single-particle Schrodinger 
wavefunctions are used to build the many-body Hamil- 
tonian with aQiD = 40. 

Below the ridge, the optimal MCBHF states are not 
energetically degenerate. The ground state and first and 
second excited states are optimized in the BHF reference 
configurations 

\N{ = 99, N 2 = 1) ground state 

\N[ = 95, N 2 = 5) first excited state (48) 

I N[ = 93, N' 2 = 7) second excited state. 



Reference configurations do not enter the Schrodinger 
based theory as the single-particle wavefunctions have 
no dependence on N[ and N 2 . It is apparent that the 
ground state is described by a binomial distribution of 
states peaked around Ni = 50 and N 2 — 50 in the basis 
where Ni and N 2 are the occupation numbers of states 
that are left- and right-localized in the trapping poten- 
tial. However, by making a unitary transformation di- 
rectly on the coefficients of the GCSs taking them back 
to the symmetric and antisymmetric basis, the expansion 
coefficients are tightly peaked around a single state. The 
optimal MCBHF ground state becomes 

1*^; {99*, 1*}) « |100, 0; {99, 1}) (49) 

while a similar result holds for the Schrodinger case. 
The underlying BHF single-particle wavefunctions as- 
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sociated with i|49|) are displayed in the first panel of Fig- 
urell3l It is important to note that these BHF wavefunc- 
tions depend strongly on the number of atoms in each 
single-particle state as well as on the interaction strength 
and do look quite different from the corresponding solu- 
tions of the single-particle Schrodinger equation in the 
second panel of Figure El Being that the MCBHF states 
(|49|l are made up of a single configuration with almost 
all N = 100 atoms in \xi), it is not surprising that GP 
theory 0, provides a good description of the ground 
state of the simple BEC. 

A harmonic oscillator like distribution of expansion co- 
efficients continues from the ground state up to and in- 
cluding the excited states that lie within the ridge at 
approximately 470 hu> in Figure ITU1 However, above the 
ridge, the distribution of GCSs making up each optimal 
MCBHF state takes on a striking new form. Macroscopic 
quantum superpositions states emerge in the spectrum. 
These states, which are pairwise degenerate, are also 
known as entangled number states and are colloquially 
called Schrodinger cats. The lower left panel in Figure 
El depicts the distribution of GCSs in the optimal 63rd, 
64th, and 65th excited MCBHF states at a barrier height 
of 3 ftw corresponding to ckqid = 40. The lower right 
panel depicts the associated distribution of Schrodinger 
Fock states. It is found that the former three states are 
optimized in the BHF reference configurations 

\N{ = 50, N!} = 50) 63rd excited state 

\N{ = 50, = 50) 64th excited state (50) 

\N[ = 50, = 50) 65th excited state. 

The underlying BHF single-particle wavefunctions for the 
63rd excited state are displayed in the second panel of 
Figure El Once again, it is important to note that these 
BHF wavefunctions depend strongly on the number of 
atoms in each single-particle state as well as on the in- 
teraction strength and do look quite different from the as- 
sociated solutions of the single-particle Schrodinger equa- 
tion displayed in the second panel of Figure El Note that 
the Schrodinger based ground and all excited states stem 
from the same Schrodinger wavefunction. The optimal 
MCBHF excited states are of the approximate form 

|*"; {50*, 50'}>„ « /r 5 |15, 85; {50, 50}) 
± 5^185, 15; {50, 50}) 

where v = 63,64,65. The coefficients f" 5 and <7g 5 are 
meant to indicate that these states are not made up of 
just two GCSs, but rather there is a distribution of GCSs 
peaked around |15, 85; {50, 50}) and |85, 15; {50, 50}) as 
is evident from the lower left panel of Figure El These 
states are, therefore, multiconfigurational in nature. In 
Figure El the "+" combination is for v = 64 and the 
"— " combination is for v = 63, 65. Sharper and more 
extreme macroscopic superposition states do appear in 
the spectrum at even higher lying energies. Notice that 
many more configurations are involved in the lower left 
panel of Figure IT21 than Fock states in the lower right. 



In addition, we point out that in the large barrier 
height limit, the BHF reference configuration underlying 
the optimal MCBHF ground state and low lying excited 
states approaches \N[ — 50, N% — 50). This can be seen 
by comparing trends in the lower two panels of Figure El 
which correspond to barrier heights of 9 htu and 13 hw. 
The distribution of GCSs at these barrier heights is also 
sharply peaked at iVi = 50 and N2 = 50. Therefore, the 
optimal MCBHF ground state takes the form 

1^; {50*, 50'}) « [50, 50; {50, 50}) (52) 

in the large barrier height limit, which is a fragmented 
BEC state. Since the optimal BHF reference configura- 
tion has approximately 50 atoms in each of two single- 
particles states, it is not surprising that a two-single- 
particle state mean field theory provides a realistic de- 
scription of the fragmented BEC ground state 0, |||. 
The MCBHF state (|52H is, after all, a single configura- 
tion with about N/2 atoms in each of two single-particle 
states. 

We have demonstrated that the optimal MCBHF 
ground states at both zero and large barrier height are 
described by a single GCS. The low barrier limit is char- 
acterized by only one single-particle state. This is why 
ground state zero temperature properties of the simple 
BEC are well understood with GP theory alone. The 
large barrier height limit requires, at minimum, two 
single-particle states. GP theory, being a one-single- 
particle state theory fails in this case, however, the two- 
single-particle state BHF theory provides a good descrip- 
tion here. Furthermore, it has been shown that states do 
exist at all barrier heights that can only be described by 
superpositions of multiple configurations into which the 
effects of atom-atom interactions are incorporated. Com- 
parison with the Schrodinger model of Section II demon- 
strates a marked difference in the atomic structure of the 
condensate due to the neglect of mean field effects on the 
underlying single-particle wavefunctions. Therefore, the 
full effects of the mean field as well as correlation must be 
included in order to gain understanding of the entire pro- 
cess of BEC fragmentation and not just an understanding 
of the simple or fragmented ground state. 

VI. CONCLUSION 

Multiconfigurational Hartree-Fock theory has been for- 
mulated for the many-body problem associated with a 
gas of identical bosonic atoms trapped at zero tempera- 
ture in potentials that can be continuously deformed from 
single to double well geometries. A didactic survey of our 
approach has been presented which clarifies many of the 
principles and approximations that are found in other 
relevant approaches from the literature. In the extreme 
limit of a single configuration, MCBHF theoryrecovers 
two-single-particle state mean field theory which 
includes the effect of the condensate mean field on the 
single-particle states but lacks the correlation that arises 
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FIG. 13: (Color online) BHF single-particle wavefunctions xi an d X2 versus coordinate z at a barrier height of 3 hui. The first 
panel corresponds to the BHF configuration that optimizes the MCBHF ground state. That is, the first panel depicts the BHF 
configuration [99, 1) associated with oqid = 40. The second panel correspond to the BHF configuration that optimizes the 
63rd MCBHF excited state. That is, the second panel depicts the BHF configuration 1 50, 50) associated with qqid = 40. All 
wavefunctions have been set at the chemical potentials /in and associated with each configuration. At each barrier height, 
we plot the corresponding Thomas-Fermi wavefunction as a solid black curve. 



between configurations. In the opposite extreme limit, 
our approach recovers the exact diagonalization of the 
many-body Hamiltonian in a restricted two-state basis of 
Fock states ^ij . Atomic correlation arises automatically 
in this case, but the full effects of atom-atom interaction 
are missing. In particular, the underlying single-particle 
states have no dependence upon the mean field of the con- 
densate. By incorporating both approaches, the MCBHF 
theory is capable of describing general many-body states 
of the system that are made up of a superposition of many 
configurations into which mean field effects are included. 

MCBHF theory has been implemented in a systematic 
study of the many-body atomic structure of the BEC 
and its fragmentation in a double well trapping potential. 
Nonlinear and linear optimization has been utilized in 
conjunction with the variational and Hylleraas-Undheim 
theorems to find the optimal ground and excited states of 
the condensate throughout its break up. A variety of in- 
teresting delocalized and localized, single configurational 
and multiconfigurational states have been found to arise 
throughout the spectrum at all barrier heights. Contri- 
butions from the condensate mean field on the underlying 
single-particle states as well as correlation effects, which 
arise between configurations, have been emphasized to 
be essential in order to obtain a more complete micro- 
scopic understanding of the condensate's atomic struc- 
ture throughout the fragmentation process. Future work 
will be devoted to the construction of an associated dy- 
namical theory that incorporates the rich atomic struc- 
ture of the MCBHF formalism. It is intended to apply 
such a dynamical theory to the accurate modeling of the 
recent experiments 0, 0, H in which the observables 
are time-dependent condensate densities, rather than sta- 
tionary state energy levels. However, the importance of 



mean field effects illustrated here for the energy levels 
of the many-body system as a function of barrier height 
clearly indicate that this same mixing of mean field and 
correlation effects will be important there. 
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APPENDIX: HYLLERAAS-UNDHEIM 
THEOREM AND OPTIMIZATION OF BOUNDS 
FOR EXCITED STATE EIGENVALUES 

A nonrigorous statement of the Hylleraas-Undheim 
theorem [51j is as follows: The N energy eigenvalues E v 
of the Hamiltonian Hki = (gk\H\gi) represented in the 
orthonormal basis {gi, . . . , gx} can be ordered so that 

Ei < E 2 < ■ ■ ■ < E N . (A.l) 

Each E v is an upper bound to the exact vth. eigenvalue 
i?° x of the same symmetry, i.e., 

E™ < Ei,E? < E 2 , ■ ■ ■ , E C N X < E N . (A.2) 
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The bounds can be systematically improved as N is in- 
creased since the energy eigenvalues necessarily move 
downward or stay the same. For example, the eigenvalues 
E' v of the Hamiltonian represented in the orthonormal 
basis {g\, . . . , g^, gN+i}, which contains one extra basis 
function, interlace the eigenvalues E v in such a way that 

El* < E[ < Ei, £f x < E 2 < Ei , ■ ■■ , Eif < E' N < En ■ 

(A.3) 
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